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Abstract 

The determination of a differential equation underlying a measured time series 
is a frequently arising task in nonlinear time series analysis. In the validation of 
a proposed model one often faces the dilemma that it is hard to decide whether 
possible discrepancies between the time series and model output are caused by an 
inappropriate model or by bad estimates of parameters in a correct type of model, or 
both. We propose a combination of parametric modelling based on Bock's multiple 
shooting algorithm and nonparametric modelling based on optimal transformations 
^ . as a strategy to test proposed models and if rejected suggest and test new ones. We 

Z!^ I exemplify this strategy on an experimental time series from a chaotic circuit where 

^^ ' we obtain an extremely accurate reconstruction of the observed attractor. 

§ ■ PACS: 05.45.-a, 05.45.Tp, 84.30.-r 
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1 Introduction 



It is an old goal in nonlinear time series analysis to infer the "Equations of 
motion from data series" ([|). Especially, for continuous flow systems modelling 
a sampled time series by a differential equation might allow for insight into the 
mechanisms at work by interpreting the resulting structure of the equation and 
values of the parameters. This is known as "interpretability" of a model (H) 
as opposed to black-box approaches like an attractor reconstruction (H). 
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The straightforward procedure to estimate the parameters in differential equa- 
tions is to estimate time derivatives from the data and determine the param- 
eters by a least-squares minimisation (^; ^; ^. This approach is firstly limited 
by the additive observational noise which usually covers the observed dynam- 
ics and prohibits the reliable estimation of the derivatives (^, especially if only 
one component of the multidimensional system can be measured. Secondly, it 
is hampered by the huge number of possible nonlinear models that have to be 
compared. 

Fortunately, there is often prior knowledge that gives constraints on the model 
or even suggests a specific type of model (|[ ||). The validity of a specific 
model can be evaluated by comparing properties of simulated time series with 
the measured one @. This approach faces a dilemma. For the simulation the 
parameters have to be specified. Thus, it is difficult to decide whether possible 
discrepancies between the properties of the simulated and the measured time 
series are caused by the fact that the chosen type or structure of model is wrong 
or by the circumstance that the parameters have been chosen inappropriately 
in the correct type of model, or both. In this paper we show that parametric 
modelling based on Bock's multiple shooting algorithm (p!0| ; |Tl|) can solve this 
dilemma. If the model type is rejected we propose a search in structure space 
instead of parameter space. A search in structure space can be performed in 
two ways: Fitting coefficients for certain basis functions (e.g. coefficients of 
polynomials) and nonparametrically. Here we chose nonparametric modelling 



based on optimal transformations and maximal correlations (JT^; |T3|; [T^) as an 
exploratory tool to suggest a new type of model that again can be tested by 
parametric modelling. We explain and exemplify this strategy on a measured 
time series from a chaotic circuit (^; |T5p. In this application our strategy yields 
a reconstruction of the observed attractor of unprecedented accuracy. 

The paper is organized as followed: In the next section we briefly describe 
the two modelling procedures. For detailed discussions of the mathematical 



methods, proofs of convergence and numerical details, see ([1^; |TT]; |1^). In 
Section 3 we introduce the data and the model derived from prior knowledge 
for these data in (^. The parametric modelling based on the suggested model, 
the search for a better model by a nonparametric procedure and the final 
parametric fit is presented in Section 4. 



2 Methods 

2.1 Bock's multiple shooting algorithm for parameter estimation 

A common setting in modelling time series by differential equations is 



f(t) = /(f(t),p) (x(t)Gi?'", peR^, te[to,tj]) (1) 

y{ti)=g{x{U))+f]iU), (2) 

where / defines the dynamics that depends on a parameter vector p. The 
state vector x{t) is observed through a function g{.) and the observation y{ti) 
is sampled at times U and disturbed by white observational noise r/(tj) of 
standard deviation a^. In general, the observation function g{.) can also contain 
unknown parameters. In the following, for ease of clarity, we assume the often 
met condition of a known, scalar observation function 

gi-) = xiit,) (3) 

which records one, for ease of notation the first, component of the dynamical 

state. 

A first approach to estimate the parameters without the need to estimate 
derivatives from the data is the initial value approach (|T^; |T^. For this proce- 
dure, initial guesses for the parameters and the initial values x{to) are chosen. 
Then Eq. (1) is solved numerically and estimates y{ti,p,x(to)) are calculated 
by Eq. (3). The error cost function 

2 _ ^ ^ {y{ti,p,x{to))-y{ti)f 
^^ i=i "* 

is minimized with respect to the parameters p and the initial values x(to) 
by some numerical optimization algorithm (|l^). For this procedure, the only 
information from the measured time series that enters the initial guesses of 
the optimization procedure is the value of the observed component, Eq. (3), 
at time to. 

Simulation studies have shown that, for many types of dynamics, this approach 
is numerically unstable by yielding a diverging trajectory or stopping in a 



local minimum ([19| ; |20| ; pi|). The reason for this is that even for slightly wrong 
parameters, the trial trajectory looses contact to the measured trajectory. This 
is most evident in the case of chaotic dynamics, where due to the sensitivity 
with respect to initial conditions the numerical trial trajectory is expected 
to follow the measured trajectory of the system only for a limited amount of 
time. This divergence of the numerical and measured trajectory introduces 
many local minima in the landscape of the error functional, Eq. (4). 

This problem can be circumvented by a multiple shooting algorithm intro- 
duced by Bock ( p!0| ; |Tl|). Here we only briefly explain this algorithm. The basic 



idea of the algorithm is to start the optimization with an only piecewise con- 
tinuous trajectory that stays close to the data. If the observation function is 



g{.) = xi{t), more information than only the first value of the measured time 
series as in the initial value approach can be used as initial guesses for the 
optimization procedure by the following strategy: The time interval [to?^/] of 
measurement is divided into numerous segments [tj,tj+i]. A trial trajectory 
for each segment is calculated using the information xi{tj) = y{tj) from the 
measured time series and initial guesses for the remaining components oix{tj). 
The condition that the underlying trajectory is smooth enters into the algo- 
rithm by a constraint in the cost function, Eq. (4). This constraint is nonlinear 
in the parameters but enters the optimization strategy using a GauB-Newton 
procedure only in a linearized way. Therefore, the trajectory is allowed to be 
discontinuous at the beginning of the optimizing iterations but is forced to 
become smooth in the end. After convergence the algorithm also provides an 
estimate x(t) of the unobserved dynamical state. 

For time series of chaotic systems it will, in general, not be possible to find 
a trajectory x{t) that shadows the true trajectory for arbitrary long times. 
Furthermore, if the model is not correct, experience shows that a fit to the 
whole time interval [to,t/] of measurement often does not converge. For both 
reasons, the time series is cut into pieces of equal length, and the parameters 
are fitted to all the pieces simultaneously. The length of these pieces is chosen 
as large as possible. 

We show the process of convergence for the time series under investigation in 
Section 4.3, Fig. 6. Analogous illustrative applications to simulated time series 



are given in (^; |22|; |2^). Note that after convergence the algorithm is identical 
to an initial value approach: It predicts the time series for each of the whole 
pieces based on the estimated initial values of the unobserved state vector. 

After convergence confidence intervals for the parameters can be calculated 
from the curvature of the cost function ([n|). Note that both approaches do 
not need an attractor reconstruction by a delay embedding. Thus, all problems 
associated with the delay reconstruction of a chaotic and noisy phase space, 
like finding an optimal embedding window and the "curse of dimensionality," 
are of no or minor importance here. Most importantly, we do not need such 
huge amounts of data as generally needed for a useful delay reconstruction 



and the method is also applicable to transient time series, see (|20| ; ^J). 



By Bock's multiple shooting algorithm the probability of stopping in local 
minima is reduced compared to the initial value approach. Nonetheless, the 
algorithm should be applied with different initial guesses for the parameters 
and the unobserved components at times tj to yield confidence in the global 
optimality of the resulting estimates. 



2.2 Nonlinear regression and optimal transformations 



Optimal transformations and the associated concept of maximal correlation 
provide a nonparametric procedure to detect and determine nonlinear rela- 
tionships in data sets. Let X and Y denote two zero-mean data sets and 

R(X,Y)- /l^^l (5) 



their (normalized) linear correlation coefficient, where E[] is the expectation 
value. The basic idea of this approach is to find transformations Q{Y) and 
$(X) such that the absolute value of the correlation coefficient between the 
transformed variables is maximized. That is, the so-called maximal correla- 
tion dl; 11; 13) 



^{x,Y)= sup |i?(e*(y),$*(x))| (6) 



e*,$* 



is calculated. The functions Q{Y) and ^{X) for which the supremum is at- 
tained are called optimal transformations. Generalizing the concept of linear 
correlation, where the linear correlation coefficient R{X, Y) quantifies linear 
dependences 

Y = aX + 7] {aeR) , (7) 



\E'(X, y) quantifies nonlinear dependences of the form 

e{Y) = ^{X) + 7] . (8) 

Especially, if there is complete statistical dependence (^), i.e., y is a function 
of X or vice versa, the maximal correlation attains unity. This is also true for 
the relation (8) with ry = 0. Here we are mainly interested in the estimation 
of the optimal transformations for the multivariate regression problem 

e(r) = $i(Xi) + . . . + ^^{x^) + v- (9) 

This is an additive model for the (not necessarily independent) input vari- 
ables Xi, . . . , Xm- The regression functions involved can be estimated as the 
optimal transformations for the multivariate problem analogous to Eq. (6). 
To estimate these in a nonparametric way, we use the Alternating Conditional 
Expectation (ACE) algorithm (0). This iterative procedure is nonparamet- 
ric because the optimal transformations are estimated by local smoothing of 



the data using kernel estimators. We use a modified algorithm^ in which the 
data are rank-ordered before the optimal transformations are estimated. This 
makes the result less sensitive to the data distribution. We remark that op- 
timal transformations for multiplicative combinations of variables Xi, . . . , X/ 
can be estimated by forming 

X, = h{X^,...,Xi), (10) 



where the hi are arbitrary functions. 

With respect to the analysis of data from nonlinear dynamical systems, the 
maximal correlation and optimal transformations have been applied to identify 
delay (inp and partial differential equations (PSf). In application to differen- 



tial equations, time derivatives have to be estimated from data. The effects 
of the noise on the estimated optimal transformations is not yet completely 
understood. On the one hand, for neglectable amounts of noise this approach 
has successfully been applied to experimental physical data of different ori- 
gin ( p9| ; pOD , yielding also quantitatively accurate results. On the other hand, 
if noise contamination is strong, this method should more be used as an ex- 
ploratory tool in the process of model selection. To minimize the influence 
of the noise on the results, the variable with the best signal-to-noise ratio, 
usually the undifferentiated time series, should be chosen as Y . In Section 4.2 
we show an application to a measured time series. 



3 The data 



The time series that we will analyze in Section 4 by the methods described 
in Section 2 was taken from an electric circuit in a chaotic regime. Technical 
details of the circuit are given in (|T3p. The data were measured at The Institute 
for Nonlinear Science of UCSD and are available in the scope of the Y2K 



Benchmarks of Predictability competition at [http : //y2k . maths . ox . ac . uh 



The model proposed in (|^; |T3p to describe the circuit reads in dimensionless 
units 



x = y 

y = —X — Sy + z (11) 

i = 7(a/(x) - z) -ay , 



A MATLAB- and a C-program for the ACE algorithm can be obtained from the 



authors or from the web page [http : //www . f dm . uni-f reiburg . de/~hv/hv . html 



where x corresponds to a voltage, y to a current and z to another voltage. The 
parameters correspond to combinations of an inductivity, the resistances and 
the two capacitances of the circuit. For the chosen experiment the proposed 
parameters are a = 15.6, 7 = 0.294, a = 1.52, and S = 0.534 (|). 



The nonlinearity is given by 



/(^) 



0.528 ifx<-1.2 

x(l-x2) if - 1.2 < x < 1.2 . (12) 

-0.528 ifx>1.2 



The measured time series corresponds to the x-component of the differential 
equation. 

The nonparametric modelling by optimal transformations requires a reformu- 
lation of Eq. (11) as a scalar higher order differential equation. The equivalent 
description reads: 

x^'^^ + {6 + 7)x + (1 + 7(5 + a)x + 7(x - af{x)) = . (13) 



For the sake of clarity, we define a = S + '~f,b = l + ^6 + o' and g{x) = 7(2; — 
af{x)). For the proposed parameters, a = 0.828 and b = 2.677. Furthermore, 
to facilitate the comparison of the results of the different approaches in Section 
4.3, based on Eq. (13), we transform the original system (11), into 



x = v 

v = w (14) 

w = —aw — bv — g{x) 



Visual inspection of the time series indicates that the variance of the obser- 
vational noise is rather small. This is confirmed by spectral analysis which 
shows a low power flattening for high frequencies (|3TD. Thus, we assume that 
the observational noise is due to discretization during sampling and follows a 
uniform distribution: 

r;~ ?7(-0.0005, 0.0005) . (15) 



4 Results 



Figure 1 shows attractor reconstructions for the measured and a simulated 
time series based on Eq. (11), and the proposed parameters given in the previ- 
ous section. For the simulated time series, the attractor is smaller, the "holes" 
in the two loops are larger and the distances between the inner edges of the 
loops and the region in phase space where the trajectories change the loops 
are smaller. 



4.1 Parametric modelling I 



To compare the measured and the simulated time series in the time domain 
we apply a restricted version of Bock's algorithm. As outlined in Section 2.1 
the algorithm finally yields estimates for the parameters as well as for the 
unobserved dynamical state (in Eq. (11) denoted by x,y,z). Here, we fix the 
parameters during the optimization process to the proposed ones and only 
allow the estimation of the state to be optimized. Figure 2 shows a segment 
of the measured time series (dotted line) and the estimated time series condi- 
tioned on the parameter values a = 15.6, 7 = 0.294, a = 1.52, and 6 = 0.534 
(broken line). We have reproduced identical results for numerous different 
choices of the initial guesses for the unobserved components y and z. The 
maximal length of the pieces applied in the estimation procedure that yielded 
a convergence was 1280/is corresponding to 64 data points, which represents 
approximately 1.5 rotations on the loops. The mean squared error between 
the measured and the estimated time series is 2.312 * 10~^. Now we are in the 
situation mentioned in the Introduction, Section 1: The proposed model type 
with proposed specific parameters is not able to reproduce the measured time 
series. 

To decide the question if a change of the parameters is sufficient to explain the 
measured data or whether a different model type has to be chosen, we apply 
the full version of Bock's algorithm. As initial guesses for the parameters we 
chose those originally proposed and also tested different ones in the same 
order of magnitude. As initial guesses for the initial values of the unobserved 
y, z-components we chose values between —10 and 10. The result is given 
in Fig. 2 (solid line). The estimated parameters are a = 1333 ± 101, 7 = 
0.00291±0.00031, a = 1.420±0.075, and S = 0.799±0.075. The mean squared 
error between the measured and the estimated time series is 7.689 * 10^^. 
The mean deviation to the measured time series is decreased, but there are 
still systematic discrepancies. This can be observed by the mismatch at the 
minima and maxima of the time series. Moreover, as in the first attempt, the 
maximal length of the pieces applied in the estimation procedure that yielded 



a convergence was 1280 /xs. As mentioned in Section 2.1 this gives evidence 
that the proposed model type is incorrect. 



4.2 Nonparametric modelling 



Figure 2 shows that there is a systematic deviation of the estimated time 
course from the measured one even for the best fit parameters. This calls for an 
alteration of the model. Since there is only one nonlinearity in the suggested 
model, see Eq. (11), we suspect that the proposed functional form of this 
nonlinearity might not be correct. To obtain a nonparametric estimate of the 
functional form, we apply the method of optimal transformations, described 
in Section 2.2. Therefore, in order to make all variables accessible, we use the 
equivalent one-dimensional third order system (13). As input variables for the 
multiple nonparametric regression problem (9) we use 

Y = x, Xi = x, X2 = X, X3 = a;(^) . (16) 



The time derivatives entering Xi, X2 and X^ have to be estimated from the 
data. The observational noise on the data is rather small, see Fig. 2 (solid 
line), presumably only resulting from the analog-to-digital conversion. Despite 
the small variance the estimated time derivatives are rather noisy. Figure 3a 
shows segments of the estimated derivatives of first to third order. The third 
derivative appears to be useless for a further analysis. However, by a proper 
filtering in frequency domain, also third-order derivatives can be recovered 
easily and with high precision, as shown in Fig. 3b. 

According to Eq. (13) we expect the optimal transformations of Xi, X2 and 
X3 each to be linear. The optimal transformation of Y should turn out to be 
a linear combination of the unknown nonlinearity in the circuit and a linear 
function. It corresponds to g{x) = 7(3; — af{x)) in the original system. 

Applying the nonparametric regression analysis, we get the optimal transfor- 
mations as displayed in Fig. 4. In estimating conditional expectation values 
as necessary in the ACE algorithm, we choose a rectangular smoothing ker- 
nel that averages over 41 data points. The value of \E'(F, Xi, X2, X3) is 0.989 
indicating that Y can be well described by the chosen variables {Xi,X2,Xs). 
The results for Xi and X3 within the 5 to 95 % quantiles of the data are 
well consistent with the expected linear behavior from Eq. (13). The optimal 
transformation for X2 can still be described fairly well by a linear function. 
Linear regression yields the parameter estimates a = 0.7963 and b = 2.5041. 
The optimal transformation for Y = x, Fig. 4(a), deviates qualitatively from 
the proposed form g{x) = 7(0; — af{x)) which exhibits a discontinuous first 
derivative at x = ±1.2, see also Fig. 7 below. It suggests a description of the 



term g{x) by a polynomial 

g{x) = Y.c,x\ (17) 



m 



i=0 



We fitted polynomials of increasing order to the nonparametric estimate of the 
nonlinearity by the optimal transformation. An order of 7 yields a fit that is 
essentially unaffected by a further increase of m, see Fig. 4a (smooth line). The 
resulting parameters are Cq = 0.1033, ci = —4.9573, C2 = —0.2462, cs = 6.887, 
C4 = 0.1744, C5 = -2.656, cq = -0.040, and cy = 0.3652. 

Figure 5 displays the reconstructed attractor based on a simulation of the 
model described by the estimated optimal transformations. The size of the 
attractor coincides with the attractor reconstructed from the measured time 
series, but the appearance differs. Since, unlike in the parametric approach, we 
have not optimized the dynamics of our model but performed only a nonlinear 
regression analysis of the variables {Y,Xi,X2,Xs), it can be expected that 
this result could still be improved. This would require, however, an extensive 
search for optimal parameters in the estimation of derivatives and conditional 
expectation values in the ACE algorithm, since the effect of noise in these 
steps is not completely understood. Rather to do that, we take this result 
as an exploratory approximation that yields useful guesses for the functional 
form of the nonlinearity and the initial parameters in parametric modelling 
again. 



4-3 Parametric modelling II 



In the following, to allow for a comparison of the different approaches, we 
use the writing of Eq. (14), where the proposed nonlinearity reads g{x) = 
lix-af{x)). 

The coefficients of the even-order monomials ffited to the optimal transfor- 
mation in the previous section are rather small. Therefore, we suspect that 
they are consistent with zero. Based on the suggested form for the nonlinearity 
from the previous section, we now apply Bock's algorithm using a sum of odd 
monomials up to seventh order as nonlinearity 



9{x) = '^C2i-lX 



2i-l 
C2i-lX 

1=1 



For the new model we were able to ffi the parameters from pieces of length 
1024 points. Figure 6 shows the ffist half of one such piece of the measured 
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time series (dotted line) and the process of convergence of Bock's multiple 
shooting algorithm (solid line). 

The convergence for rather long pieces as well as the coincidence of the finally 
estimated time course of x(tj) and the measured data in Fig. 6c suggests that 
the polynomial nonlinearity of Eq. (18) provides a better description of the 
underlying system than the saturating nonlinearity of Eq. (12). 

As mentioned in Section 3, the small observation noise on the data leads to 
unrealistic small confidence intervals for the estimated parameters. To obtain 
a realistic impression of the variability of the parameter estimates, we divided 
the time series in 8 parts and report the mean and standard deviations cal- 
culated from the results for these parts: a = 0.759 ± 0.021, b = 2.526 ± 0.014, 
ci = -4.56 ±0.16, C3 = 6.44 ±0.58, Cg = -2.55 ±0.41, and C7 = 0.366 ±0.074. 
The mean squared error is 1.581 * 10^^. 

Figure 7 compares the proposed nonlinearity g{x) = j{x — af{x)), the optimal 
transformation Q{Y = x) and the corresponding function based on the final 
parametric fit. Interestingly, the fitted polynomial closely follows the proposed 
nonlinearity for small absolute values of x, but approaches the estimated op- 
timal transformation for larger values of x. 

Analogous to Fig. 1, Fig. 8 displays the reconstructed attr actor based on a 
simulated time series using the polynomial Eq. (18) as nonlinearity and the 
best fit parameters. The attractors are in excellent agreement. 



5 Discussion 



We proposed a three-step procedure for modelling nonlinear time series by 
differential equations and exemplified this strategy on a physical application. 
We chose the algorithmically more difficult modelling by differential equations 
in favor to discrete-time difference equations because the results of the former 
are usually easier interpretable in terms of the underlying physics (pBf). 



Bock's multiple shooting approach to parameter estimation in differential 
equations does not require to estimate time derivatives from the data which 
limits the applicability of many other approaches. Note that even for the very 
clean data in our application it is not possible to estimate the third deriva- 
tive without using some kind of low-pass filtering, see Fig. 3. The price to be 
payed in Bock's algorithm is that this approach needs a parameterized model. 
On the other hand, nonparametric modelling by optimal transformations does 
not require a specific parameterized model but is more susceptible to noise. 
Furthermore, it treats the problem as a case of regression, not taking into 
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account that the data were generated by a dynamical system. Nevertheless, 
as our application has shown, it can be applied to inspire parametric models 
that, again, can be checked by parametric modelling of the dynamical system 
by Bock's multiple shooting algorithm. 

The difference between the results for the nonparametric and the final para- 
metric analysis displayed in Fig. 7 shows that a simple parametric fit to the 
nonparametric estimate of the nonlinearity as reported in Section 4.2 can be 
improved by the third step of our procedure that evaluates the predictive 
ability of the model in the time domain. In the discussed application our final 
result is extremely accurate considering that we used an experimental time 
series and that we can reproduce the global dynamics of this highly nonlinear 
system (Fig. 8). The excellent agreement between the dynamics of a measured 
chaotic time series and an estimated model is a highly nontrivial result ( P2| ) 
in nonlinear modelling. 

Canonically, dynamical systems are given as vector-valued first order differ- 
ential equations. A precondition to apply the second, nonparametric, model 
generating step of the proposed strategy is that the dynamical system can be 
expressed as a scalar, higher order differential equation. It follows from The- 
orem III of Takens' seminal paper ( |5BD that for an m-dimensional system of 
first order differential equations there is always an equivalent one dimensional 
differential equation of maximum order 2m + 1. Unfortunately, it depends on 
the given system whether it is possible to find an explicit form for the one 
dimensional counterpart. Fortunately, as far as known to the authors, for all 
chaotic standard systems the one-dimensional writing is possible. Note that 
in the investigated case of electronic circuits (due to Kirchhoff's laws (^4[)) a 
huge class of nonlinear circuits can be modelled by differential equations like 
Eq. (11) anyway. 

It has to be emphasized that the suggested procedure of alternating paramet- 
ric modelling by Bock's algorithm and nonparametric modelling by optimal 
transformations is not a general purpose procedure in the sense of "Equations 
of motion from a data series" (|l|). The success of the discussed application 
depended on prior knowledge about the system, i.e. a roughly correct first 
suggestion on the structure of right hand side of the differential equation. 
This, however, is often the case in systems where the dynamics is qualita- 
tively understood, and one is interested in obtaining coefficients or special 
forms of nonlinearities involved in the dynamics. Therefore, we assume that 
our approach will be applicable mainly for systems from physics, like electronic 
circuits or lasers (§ ^T]), engineering, e.g. effects of nonlinear friction (|35|) , bio- 
chemistry, e.g. dynamics of protein folding (|36D, and biophysics, e.g. dynamics 
of photosynthesis (^4]). 
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Fig. 1. Reconstructed attractors. (a) Prom the measured time series, (b) From 
the simulated time series based on Eq. (11) and parameters a = 15.6, 7 = 0.294, 
a = 1.52, and S = 0.534. The delay time is 15 sampling units, corresponding to 
300/xs. 




Fig. 2. Segment of the measured time series (dotted line), the best fit trajectory with 
fixed parameters a = 15.6, 7 = 0.294, a = 1.52, and 5 = 0.534 (broken line) and 
result for the best fit parameters a = 1333, 7 = 0.00291, a = 1.420, and S = 0.799 
(solid line). The fit trajectories contain eight continuous pieces of 64 data points 
each. 
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Fig. 3. Piece of the time series and, from top to bottom, pieces of first to third 
time derivatives estimated from the time series. To obtain the highest accuracy, 
the derivative estimations are performed in frequency domain by multiplying the 
transformed data with the first to third power of the wave numbers, respectively, 
and subsequent back-transforming into the time domain, (a) is a direct estimate, 
(b) uses a cut-off at 10% of the Nyquist frequency corresponding to a low-pass filter. 
The reliability of the derivative estimates can be checked roughly by noting that 
the extrema of the ith derivative (i = 0,1,2) always correspond to a zero in the 
{i + l)th derivatives one line below. 
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Fig. 4. Estimated optimal transformations for the one-dimensional third order 



differential equation Eq. (13). (a) 6(1" 
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(d) $3(^3 



"(3)1 



(b) <5i(Xi 
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In (a) also a fit of a 7th order polyno- 



mial is shown. The coefficients of the nonlinearity g{x) estimated this way are given 
in the text. The vertical lines indicate upper and lower 5 % quantiles of the data. 
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Fig. 5. Reconstructed attractor based on the results for the nonparametric modelling 
by the optimal transformations. 
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Fig. 6. Segment of the measured time series (dotted line) and convergence of the mul- 
tiple shooting approach based on the polynomial nonlinearity (solid line), (a) Trial 
trajectory for initial guesses of the parameters and initial values, see text for details, 
(b) Trial trajectory after 3 iterations, (c) Result after convergence. 
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Fig. 7. Comparison of the proposed piecewise differentiable model nonlinearity, 
g{x) = 7(x — af{x)) (bold line), the nonparametrically given optimal transfor- 
mation Q{x) (dithered line), and the parametrically fitted odd polynomial (dotted 
line). 





Fig. 8. Reconstructed attractors. (a) From the measured time series, (b) From the 
simulated time series based on Eqs. (14,18) and the best fit parameters, see text. 
The delay time is 15 sampling units, corresponding to 300^s. 
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